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I. INTRODUCTION 

In recent years the nonequilibrium dynamics of low- 
dimensional quantum systems has been increasingly at- 
tracting the attention of experimentalists and theoreti- 
cians from different areas of physics An excit- 
ing field in which great experimental progress has been 
achieved within the last decade is the one of ultracold 
quantum gases. There, advances in atom waveguide 
technology 0, H, [1, [f| Q, the realization of quantum 
gases in very anisotropic traps 0, @, H|, and loading 
Bose-Einstein condensates (BEC's) on optical lattices 
[M GH, 03, GS 0, EE M, [l3 have allowed experimen- 
talists to obtain a wide variety of systems in which the 
effects of the reduced dimensionality can be studied in 
detail. 

Due to the very high control than can be achieved in 
an experiment with a degenerate quantum gas, one can 
prepare it under very specific initial conditions and study 
its evolution. This has been done in the one-dimensional 
(ID) regime in various experiments. For example, the 
transport properties of ID Bose gases on a lattice have 
been studied in Refs. [H|,[l(|, where the gas was displaced 
from the center of the trap and allowed to oscillate. More 
recently, Kinoshita et al. fl7| have addressed experimen- 
tally the question of whether an isolated integrable or 
nearly integrable system can relax to the thermal equi- 
librium state [HI]. 

On the theoretical side, integrability in low- 
dimensional systems allows one to perform exact studies 
of the equilibrium properties and the nonequilibrium dy- 
namics of well-known models pjl, H, Hi], H, HI| , some of 
which have already become relevant to experiments. One 
particular model in which we are interested in this work 



is the one of impenetrable bosons in ID [Tj, LLa, 53 ■ ^ 
has been shown theoretically 0, HH, HH that in certain 
ID regimes of low densities and low temperatures, bosons 
behave as a gas of impenetrable particles also known as 
Tonks-Girardcau bosons [hard-core bosons (HCB's)]. 

The ID homogeneous gas of HCB's was introduced by 
Girardeau [27|, who also established a one-to-one cor- 
respondence (Bosc-Fcrmi mapping) between ID HCB's 
and spinlcss fcrmions. This mapping has been rece ntly 

used to study the nonequilibrium case [II, 11, [1, si a 

where the density dynamics revealed dark-soliton struc- 
tures 128||, breakdown of the time-dependent mean- field 



thco 




| , interference patterns of the thermal gas on a 



ring 3C | , and other interesting effects during the expan- 
31,13. 



Hard-core bosons have been also realized experimen- 
tally in the presence of a lattice along the ID tubes [H |. 
In the periodic case, the HCB lattice Hamiltonian can 
be map ped onto the ID XY model of Lieb, Schulz, and 
Mattis [33]. The ID XY model has been extensively 
studied in the literature, and the asymptotic behavior of 
the correlation functions is known [33, SH, 11, 13 • With 
an additional confining potential, the case relevant to 
the experiments in Ref . [lj] , this model has been studied 
be means of an exact numerical approach in Ref. [38| . 
There it was shown that one-particle correlations exhibit 
a universal power-law decay. The generalization of the 
approach in Ref. [H[ to the nonequilibrium dynamics 
[33 . |40L l4lj revealed that during the expansion of the 
HCB gas on the lattice two very different regimes can be 
identified. If the expansion starts from a Mott insulat- 
ing state, quasi-long range correlations develop between 
initially uncorrelated particles, producin g th e emergence 
of quasicondensates at finite momentum [39l . [4l| . On the 
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other hand, for low initial densities (superfluid state), 
the momentum distribution of expanding HCB's rapidly 
approaches that of noninteracting fermions [13, El| . 



In this work we study the noncquilibrium dynamics 
of hard-core bosons on ID superlattices. The superlat- 
tice is obtained adding an extra periodic potential to the 
already existing lattice. In the soft-core regime these sys- 
tems have been already realized experimentally [42L [43| 
and studied theoretically in ID by various mean-field 
and pcrturbative approaches 44j , quantum Monte Carlo 
simulations (4f|, and exact diagonalization [46[. In the 
hard-core regime, due to the mapping to noninteracting 
fermions, one can realize that the effect of the extra pe- 
riodic potential is to open gaps at the boundaries of the 
reduced Brillouin zone. This means that in addition to 
the insulating phase with density 1 (full filling) new in- 
sulating phases appear at fractional fillings. 



In three-dimensional systems the study of the half- 
filled hard-core boson model allowed one to prove rig- 
orously the existence of Bose-Einstein condensation and 
Mott insulatin g ph ases tuning the strength of the ad- 
ditional lattice [43] ■ Here we quench the strength of the 
superlattice potential to study the dynamics of these sys- 
tems in the superfluid and insulating regimes. We are 
interested in understanding how the system approaches 
equilibrium, if it does, and in testing the prediction power 
of a generalized Gibbs ensemble theory recently proposed 
in Ref. [Hj]. In order to do so, we also study the case 
in which an additional harmonic confining potential is 
introduced in the system, as relevant to most of the ex- 
perimental situations. At low densities we obtain results 
similar to that of the recent experiment by Kinoshita et 
al. [H. 



II. HARD-CORE BOSONS IN PERIODIC 
SUPERLATTICES 

The hard-core boson Hamiltonian in the presence of a 
superlattice potential can be written as 



H = - 



2ni 

ti,i+ib}k+i + H.c. ) + A ^ cos ( — 



(1) 

where the HCB creation and annihilation operators at 
site i are denoted by b\ and b i , respectively, and the local 
density operator by hi = SJS,-. In different sites the HCB 
creation and annihilation operators commute as usual for 
bosons: 

[b i ,b]} = [b i ,b j ] = [blb]} = 0, for ijtj. (2) 

However, on the same site these operators satisfy anti- 
commutation relations typical for fermions: 



1. 



St 2 = S? = 0. 



(3) 



These constraints avoid double or higher occupancy of 
the lattice sites. In Eq. (JT|), the hopping parameters are 
denoted by £i,i+i = t. The last term represents the su- 
perlattice potential with strength A and unit cells with 
L sites. 

Using the Jordan- Wigner transformation [48| 
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f)=\ 



k H'"' ./ (4) 



one can map the HCB Hamiltonian onto the one of non- 
interacting spinless fermions, 



Hp 



2m 
I — I 



',(5) 



where f! and fi are the creation and annihilation oper- 
ators for spinless fermions at site i and hj = fjf i is the 
local particle number operator. Considering that 



The exposition is organized as follows. In Sec. II we de- 
scribe the ground-state properties of HCB's in a periodic 
superlattice, paying special attention to the behavior of 
the one-particle correlations in the superfluid and insu- 
lating regimes. Section III is devoted to the study of 
the noncquilibrium dynamics of the system after a sud- 
den switch-on and -off of the superlattice potential — i.e., 
when going from the superfluid to the insulating phases 
and vice versa. Relaxation and the collapse and revival 
of the zero-momentum peak of the momentum distribu- 
tion function are studied in detail. In Sec. IV we analyze 
the ground-state properties of the system in the presence 
of a trap, which produces a coexistence of superfluid and 
insulating phases. The noncquilibrium dynamics in the 
presence of a trap is studied in Sec. V. Finally, the con- 
clusions are presented in Sec. VI. 



N 



bN b i = -f N fl ex P I ~ i7T ^2n f I . (6) 

0=1 



where N is the number of lattice sites, one notices that 
t' N 1 = exp [— m(Nf) + 1)] tjv,i; i-e., for periodic chains 
when the number of particles in the system (JVj, = 

2~2i(fii) = 2~2i ("■{)) i s °dd, the equivalent fermionic 
Hamiltonian satisfies the same boundary conditions than 
the HCB one; otherwise, if iVf, is even, the oppo- 
site boundary conditions are required for the fermionic 
Hamiltonian. 

With the help of the mapping above one can realize 
that for HCB's in a superlattice insulating phases appear 
at fractional fillings (n^ = i/L, with i = 1, . . . ,L — 1, un- 
less a band crossing occurs |45j), in addition to the full 
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filling insulator (n, = 1) already present in the absence 
of the superlattice. This is because in the equivalent sys- 
tem of noninteracting fermions the superlattice potential 
opens gaps at the boundaries of the reduced Brillouin 
zones. In what follows, for simplicity, we restrict our 
study to superlattices with L = 2. In this case the insu- 
lating phases occur at half and full filling. The general- 
ization of our results to larger values of L is straightfor- 
ward, and does not (qualitatively) change the results we 
discuss for L = 2. 

In order to study the bosonic one-particle correlations 



Pij 



= 6 6, 



(7) 



and related quantities like the momentum distribution 
function 



N ^ 



(8) 



we follow the exact approach described in detail in Ref. 
[HI . This approach allows us to study very large system 
sizes (thousands of lattice sites) in a very efficient way. 
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FIG. 1: (Color online) Decay of one-particle correlations (av- 
eraged per unit cell as described in the text) in the superfluid 
phase. These systems have up to 1900 lattice sites and L = 2. 
The plot for p — 0.25, A — 5t was displaced down for clarity. 
Straight solid lines correspond to power laws 1/y/x. The inset 
shows Uk=o (top) and n^ = ± 1I i a (bottom) vs AT(, for systems at 
quarter filling and A — 0.5t. The straight lines signal y/Nb 
behavior. Distances are normalized by the lattice constant a. 

In Fig. [1] we show how the one-particle correlations 
decay in the presence of a superlattice potential in the 
superfluid phase. In the figure we have averaged the cor- 
relations measured from the even and odd sites in order 
to minimize the effects of the different density in the sites; 
i.e., we have plotted 

Px = ^ (Pi odd .i° dd +x/a T Pi" 1 '"" ,i" v '' n +x/ci) • 

As can be seen in Fig. [TJ one-particle correlations in the 
superlattice decay with the same power law p m ~ 1/ \J~x 
that was shown to be universal in absence of the superlat- 
tice potential [3S| . The only effect that the superlattice 



introduces and that is evident only for large values of A 
is an oscillatory behavior in p x on top of the 1 / yfx decay. 

The existence of quasi-long-range one-particle correla- 
tions implies that, like in the usual lattice, a sharp peak 
in the momentum distribution function at k = signals 
the superfluid state. The additional oscillatory behavior 
seen in p x (Fig. [T]) is reflected by additional peaks in 
at ka — ±7r, which deplete the one in k = from its value 
for ^4 = 0. This can be seen in Figs.^a) and^c), where 
we have plotted the momentum distribution function for 
two different values of the strength of the superlattice 
potential. A simple calculation also allows to extract 
from p x ~ 1/y/x the scaling behavior of the k = and 
k = ±7r/a peaks as the system size is increased keeping 
the density constant. One finds that n^ = Q t ± w / a ~ y/Nb. 
Such scaling can be seen in the inset of Fig. [T] where we 
have plotted these quantities for systems at quarter filling 
and A = 0.5t. 
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FIG. 2: (Color online) Momentum profiles for periodic sys- 
tems with 900 lattice sites in a superlattice potential with 
L = 2. In (a) and (c) the system is in a superfluid state. 
The peaks at k = and k = ±rr/a signal the presence of 
quasi-long range one-particle correlations (Fig. [TJ. On the 
other hand, in (b) and (d) the system is in an insulating state 
where one-particle correlations decay exponentially (Fig. [3} . 



At half filling the behavior of the system is completely 
different to the one above. As mentioned before a gap 
opens in the spectrum. Since Eq. ^ can be easily diag- 
onalizcd for L = 2, one immediately obtains that the gap 
is A = 2A, and the dispersion relation in the two bands 
is 



e±(fc) = ±y/U 2 cos 2 (ka)+A 2 , 



(9) 



where by 
lower one. 



we mean the upper band and by "— " the 



The presence of this gap produces an exponential decay 
of the one-particle correlations as shown in Fig. [3J The 
correlation length £ is a function of the gap — i.e., of A. 
We have calculated the correlation length £ as the second 
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moment of the onc-particlc density matrix, 



1 



1 2y ( x i x j) Pij 

2 Eij Pij 



(10) 



which for large values of £ is equivalent to calculating it 
fitting an exponential decay p x ~ exp(— as shown 
at finite temperatures in Ref. [49[ . 
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FIG. 3: (Color online) Decay of one-particle correlations (av- 
eraged per unit cell as described in the text) in the insulating 
(half-filled) case. The systems considered have up to 1500 lat- 
tice sites and L — 2. The straight lines signal an exponential 
decay with a correlation length £, which is a function of A. 
In the inset we show £ vs A. The straight lines depict the 
asymptotic behavior of £. For very small values of A one has 
that £/a ~ l/(A/t) while for large values of A one obtains 
that £/a ~ 1/y/Aji (where a is the lattice constant). 

The behavior of £ as a function of A is shown in the 
inset of Fig. [3l There one can see that it exhibits two 
different functional forms for small and large values of 
A. For small values of A we find that £/a ~ l/(A/t), 
while for large values of A it is £/a ~ 1/ \J ' A/t (a is the 
lattice constant). In general, a correlation length £ ~ 1/A 
develops from free bosonic excitations across a gap ~ A. 
Such an argument results already at a mean-field level, 
and in the absence of diverging correlation lengths, it is 
qualitatively correct for small values of A. For very large 
values of A one can use perturbation theory to determine 
the asymptotic behavior of £. For A/t — > oo the ground 
state at half filling is just an array of Fock states with 
one particle in the odd sites and no particle in the even 
sites I^A/t-oo = Ili WiodJOk^n- Up to first order in 
perturbation theory the ground state of the system can 
be written as 

|*)g « l*>A/ t ^oo + ^4 2 + H ' c -) l*>A/t^oo, 

.(11) 

which means that for i ^ j [the only ones that contribute 
to the numerator in Eq. (|10p] 



t 

2A 




In the denominator in Eq. (jlpp only the terms i = j (the 
densities) contribute to first order. Hence, 

^ = a \/I (13) 

for very large values of A/t. Indeed, £/a = y/t/A is the 
straight line that in the inset in Fig. [3] follows the data 
points for large values of A. 

The consequence of the exponential decay of the one- 
particle correlations in the momentum distribution func- 
tion can be seen in Figs. El^b) and[2jd). There we have 
plotted rik for two different values of A and compared 
it with A = 0. The effect of the gap at half filling is 
dramatic. It destroys the peaks in the momentum dis- 
tribution function. This feature can be used experimen- 
tally to detect the presence of an insulating state in a 
superlattice, like it has been done previously to detect 
the presence of a Mott insulator in the absence of the 
superlattice. 

An overall picture of the behavior of the momentum 
distribution function when changing the density can be 
obtained plotting rik=o vs P as shown in Fig. [4] For con- 
trast we have also plotted the result for A = 0. Figured] 
shows that while for low densities the effect of A is almost 
imperceptible, it becomes very large as the density ap- 
proaches half filling. In the insulating state nk=o attains 
its minimum value. 




for = 1, 

for |i-j|>l. (12) 



FIG. 4: (Color online) Normalized occupation of the k = 
momentum state as a function of the density in periodic 
systems with iV = 1000. For A ^ one can see that rifc = o is 
strongly suppressed around n = 0.5, where the system is an 
insulator. 



III. DYNAMICS AND RELAXATION IN THE 
SUPERFLUID AND INSULATING REGIMES 

We study in this section the dynamics of hard-core 
bosons on superlattices when by a sudden change of the 
strength of the superlattice potential one goes from a su- 
pcrfluid to an insulating regime, and vice versa. Since for 
L = 2 this can only occur at half filling, we will restrict 
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our analysis to that case. In addition, for a closer con- 
nection with the experiments instead of periodic systems 
we consider open boxes, where translational invariancc 
is broken, but keeping the same phase diagram than in 
the periodic case. We will consider systems confined in 
harmonic traps, where superfluid and insulating phases 
can coexist, in the next two sections. 

We start our study with the case in which the initial 
state is superfluid, at half filling with A = 0, and A 
is suddenly changed to a finite value, for which in the 
ground state the system would be an insulator. This is 
similar, in the Hubbard model, to a change of the on-site 
repulsion U from a value in which the system is superfluid 
to U > U c (U c being the critical value for the formation 
of a Mott insulator). The number of particles per site in 
this case has to be integer. Such study for a Hubbard like 
experimental system has been done in three-dimensional 
optical lattices in Ref. [5011. 

In the experiment [5C| it was found that as the system 
evolves a collapse and revival of the initial momentum 
distribution function (the interference pattern) occurs. 
This can be easily understood considering that for time 
scales much smaller than the one set by the hopping pa- 
rameter, the particle number per lattice site remains al- 
most unchanged and only the phases evolve. This evolu- 
tion is dictated by the on-site interaction between atoms 
(U), which in the Hubbard language reads 

Hint = \un{n - 1). (14) 

Hence, the evolution operator has the form exp[— iUn(n— 
1)t/2], where r is a time variable (which we will give 
in units of h in what follows). This means that phases 
in different lattices sites evolve differently, according to 
the number of atoms present, collapsing the interference 
pattern. However, for times T rev = 2irn/U, where n is an 
integer, a revival of the interference pattern occurs. 

In Fig. [5] we show the time evolution of the occupa- 
tion of the zero-momentum state for a box that initially 
has A = and the evolution is performed by a sudden 
change to three different values of the final A. This figure 
shows that also in a ID superlattice loaded with HCB's 
a collapse and revival of the zero-momentum peak oc- 
curs. Here, double or higher occupancy is forbidden by 
the hard-core constraint (U = oo) so that the short-time 
evolution is only determined by the term in the Hamil- 
tonian proportional to the superlattice potential, 



H 



A)% 



(15) 



which means that A is the parameter that controls the 
revival time (equivalent to U in the Hubbard case) . This 
can be better seen in the insets in Fig. where we have 
plotted the first three revivals of rik=o- Comparing the 
times with the values of A one can see that indeed T rev ~ 
irn/A. 

The second effect that is apparent in Fig. [3] is 
the damping of the collapse and revival of the zero- 
momentum peak of rife. It is related to the change of 




FIG. 5: (Color online) Evolution of the occupation of the zero- 
momentum state after A is changed from A — in the initial 
(superfluid) state to A = l.Ot (a), A = AM (b), and A = 8M 
(c). The insets show in more detail the first three revivals of 
nfc=o- Dashed lines show rik=o obtained from the constrained 
thermodynamics theory (see text). These systems have 300 
lattice sites, 150 particles, and L = 2. 



the occupation of each lattice site due to the finite value 
of the hopping parameter t. Variations in the density 
destroy the periodic evolution set by Eq. (fl~5"|) . In these 
systems the dynamics of the density is determined by a 
combination of the times scales given by t and A. (Notice 
that A in the even sites acts like a barrier between odd 
sites.) In Fig. [3] one can see that the largest damping 
occurs for the smallest the value of A/t. For the largest 
value of A shown in the figure (A = 8t) the collapse and 
revival is almost completely damped after r = lOt. 

After observing the damping of the collapse and re- 
vival of the momentum distribution function in Fig. 
one immediate question one could ask is to what kind of 
momentum distribution is the system relaxing. In nonin- 
tegrable systems one expects to relax to the thermal 
distribution. In the grand canonical ensemble it can be 
obtained using the many-body density matrix 



p=Z- L cxp - [H - fxN b ) /k B T , 



(16) 



6 



where 



1.5 



Z = Tr 



{exp [- (H - AtiV b ) /fc B r] } 



(17) 



is the partition function. In order to determine the rele- 
vant temperature (kgT) and the chemical potential (fx) 
one can use the energy and number of particles of the 
evolving system, which do not change during the dy- 
namics; i.e., kgT and \x can be calculated using the con- 
straints 



E = Tr Hp , N b = Tr N b p 



(18) 



In Fig. [S] we compare the momentum distribution func- 
tion obtained within the grand-canonical approach de- 
scribed above (called "thermal" in the figure) with the 
time-averaged momentum distribution obtained during 
the evolution, after the damping of the oscillations shown 
in Fig. [5] We have averaged the time evolving n,fc(r) be- 
cause, as can be seen in Fig. [5l this quantity exhibits 
fluctuations around its mean value. The exact njt's for 
the thermal equilibrium were obtained using the method 
detailed in Ref . [i^] . Figure[B]clearly shows that the mean 
values around which nk(r) fluctuates are not the ones of 
the usual thermal equilibrium; i.e., these systems are not 
ergodic in the usual sense. Since hard-core bosons in ID 
lattices are integrable (33j , nonergodic behavior may have 
been expected due to the presence of extra constants of 
the motion (whose number is infinite in the thermody- 
namic limit) in addition to the energy and the number of 
particles. 

The problem of relaxation in an integrable system, like 
the one of hard-core bosons that we are analyzing here, 
has been already discussed in Ref. [l8| . There it was 
conjectured that the correct many-body density matrix 
that describes the properties of an integrable system af- 
ter relaxation is given by a generalization of the Gibbs 
ensemble 



Pc 



z: 



where 



Z r . = Tr 



exp 



exp 



£a, 




(19) 



(20) 



is the corresponding partition function, {I m } is a full set 
of integrals of motion, { A m } are the Lagrange multipliers, 
and m = 1, . . . , TV. These Lagrange multipliers can be 
calculated using the expectation values of the full set of 
integrals of motion of the evolving system, which do not 
change with time 



(In 



Tr L 



(21) 



In Eq. (|2ip and in what follows, (■ • • ) T means expecta- 
tion values in the time-evolving system when they do not 
depend on time. 
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FIG. 6: (Color online) Time average of the momentum dis- 
tribution function after the damping of the oscillations seen 
in Fig. [5] We averaged measurements done in time intervals 
At = 20t between times r = 20t and r = 5000t. The time 
average is compared with the results obtained in the usual 
thermal ensemble and the constrained theory explained in the 
text. These systems have 300 lattice sites, 150 particles, and 
L — 2. In all cases the initial A = and the final one A = l.Ot 
(a), A = 4.0t (b), and A = 8.0t (c). The corresponding tem- 
peratures in the grand-canonical ensemble are k B T = 0.86t 
(a), k B T = 6.86i (b), and k B T = 25.75t (c). 



The conjecture above is still based on the ergodic 
hypothesis, but generalized to consider that the region 
available of phase space is determined by all constants 
of the motion, and not only by E and N as usual for 
nonintegrablc systems. Hence, what the density matrix 
defined by Eq. (|19p does is to maximize the many body 
Gibbs entropy, 



S = feglr 



Pc In I -!- 

Pc 



(22) 



subject to the constraints imposed by all the integrals of 
motion. 

Like in the ground state [38[ and the thermal equi- 
librium case [49|], in order to calculate the HCB expec- 
tation values using Eq. (fT9|) we take advantage of the 
Jordan- Wigncr transformation and the mapping to non- 
interacting fermions. In fcrmionic language the integrals 
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of motion can be easily constructed after diagonalizing 
the Hamiltonian ([5]): 



Hf7%\0} = E m j£\0). 



(23) 



Since the occupation of the eigenstates of the final 
fermionic Hamiltonian cannot change in time, these 
fermions are noninteracting; a complete set of integrals 
of motion can be constructed to be 



1 m I V ' m 'm J 



(24) 



The constraints in Eq. (|20p . in the fermionic represen- 
tation, lead to an analytical expression for the Lagrange 
multipliers, 



A m = In 



1 - (lL)r 



(I 



f< 



(25) 



which allows us to build the equivalent fermionic Hamil- 
tonian and to obtain the HCB expectation values us- 
ing the approach explained in detail in Ref. (49|. One 
should notice at this point that the constraints defined 
by Eq. (|24p , when written in the bosonic language, in- 
volve many bosonic creation and annihilation operators. 
Hence, these constraints lose the bilinear character they 
have in the fermionic representation. 

In Fig. [5] one can see that the results obtained with the 
constrained thermodynamics theory explained above are 
indistinguishable from the ones of the time average af- 
ter damping. Hence, even for the evolution of these pure 
quantum states in intcgrablc systems one can define a 
constrained ensemble able to predict the mean values of 
observables after relaxation. We have also included in 
Fig. [5] the value of rik=o obtained from the constrained 
thermodynamic theory for a comparison with the time 
evolving nk=o (t) . There one can see that after relaxation 
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FIG. 7: (Color online) Evolution of the occupation of the zero- 
momentum state after A is changed from A — 1.04, A = 4.04, 
and A = 8.0t, (from top to bottom) in the initial (insulating) 
state to A — 0. Dashed lines show nk=o obtained from the 
constrained thermodynamics theory (see text). These sys- 
tems have 300 lattice sites, 150 particles, and L = 2. 



the small fluctuations of nk=o{ T ) indeed occur around the 
value obtained with the generalized Gibbs ensemble. No- 
tice that the case described so far in this section seems 
to be one of the worst scenarios for HCB's since the evo- 
lution is performed with a Hamiltonian for which the 
ground state of the system is gapped; i.e., relaxation is 
expected to occur slowly due to the presence of the gap. 

In what follows we consider the inverse case. The case 
in which initially the system is in an insulating state and 
by turning off the supcrlatticc potential this state is al- 
lowed to evolve with a Hamiltonian whose ground state 
is supcrfluid. 

In Fig. [71 we show the time evolution of nk=o from 
three initial insulating states. They have the same values 
of A used during the evolution of the systems in Figs. [5] 
and [5] The evolution is performed after turning off the 
supcrlatticc potential (making A = 0). In contrast to 
the case in which the superlattice is turned on, one can 
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FIG. 8: (Color online) Time average of the momentum distri- 
bution function measured in time intervals At = 204 between 
times r = 204 and r = 50004. The time average is compared 
with the results obtained in the usual thermal ensemble and 
the constrained theory explained in the text. These systems 
have 300 lattice sites, 150 particles, and L — 2. Initially 
A = l.Ot (a), A = 4.04 (b), and A = 8.0t (c), and in all cases 
the final value of A is A = 0. The corresponding tempera- 
tures in the grand-canonical ensemble are ksT = 0.634 (a), 
k B T = 2.064 (b), and k B T = 4.034 (c). 
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see in Fig. [7] that no collapse and revivals occur in rik- 
After an increase of the occupation of n k =o the system 
steadily relaxes to an approximately constant momentum 
distribution. 

As before, in Fig. [S] we compare the results of the time 
average of the momentum distribution function after the 
system has relaxed to a stationary distribution with the 
expectation values of this quantity in the grand-canonical 
( "thermal" in the figure) and generalized ( "constrained" 
in the figure) Gibbs ensemble. Remarkably, when one 
starts from an insulating state, where correlations decay 
exponentially — i.e., there is no quasi-long-range order — 
the results obtained within a grand-canonical description 
are very similar to the ones obtained during the time 
evolution. As can be seen comparing Figs. [8fa) , EJb), 
and[5Ic), the deeper the initial state is in the insulating 
regime the better is the agreement between the grand- 
canonical and the average of the time evolving n^s after. 
Actually, only in Fig.[5]Ja) can one clearly distinguish the 
differences in the occupation of the momentum states 
with very low values of k. On the other hand, the results 
obtained within the constrained theory are in all cases 
indistinguishable from the ones obtained during the time 
evolution. 
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FIG. 9: (Color online) Comparison between the time average 
of the momentum distribution function of two systems start- 
ing from different initial conditions but having the same en- 
ergy during the evolution. (As in the previous figures, nk for 
the time average was measured in intervals At = 20t between 
times r = 20t and r = 50004.) In both cases the evolution is 
done with a Hamiltonian in which A = l.Oi, but the initial 
state is in one case a superfluid state (A — 0), and in the other 
an insulating state with A = 6.6i. The time average is com- 
pared with the results obtained in the usual thermal ensemble 
(which depends only on the final energy and the number of 
particles — i.e., is the same in both cases) and the constrained 
theory explained in the text. These systems have 300 lattice 
sites, 150 particles, and L = 2. In the legend "(TA)" means 
time average in the out-of-equilibrium system and "(C)" the 
result of the constrained thermodynamics. 

A further proof of the "memory" that these intcgrablc 
systems have on the initial conditions can be obtained 
comparing the final momentum distribution function to 



which the system relaxes starting from two different ini- 
tial states that have the same final energy. This compar- 
ison is done in Fig. [5] where we show results for two sys- 
tems that start their evolution from a supperfluid (zero 
A) and insulating (nonzero A) states, evolving with a 
Hamiltonian in which the ground state is an insulator 
(nonzero A). While the usual grand-canonical descrip- 
tion anticipates that they both should relax to the same 
momentum distribution, it can be seen in the figure that 
this is not the case. Only the generalized Gibbs distri- 
bution introduced before predicts two final momentum 
distribution functions and hence is able to describe the 
mean values of n k after relaxation. 

To conclude this section we would like to make some 
remarks about the one-particle density matrix. While 
in the usual and generalized Gibbs distributions one can 
easily work in a basis where this quantity is real, this 
is not the case during the noncquilibrium dynamics of a 
quantum system. In the latter case, in all our calcula- 
tions, the one-particle density matrix is a complex ob- 
ject (Px{t) = \Px(t)\ exp[— i9 x (T)]) and nontrivial time- 
dependent phases [#;c(t)] enter into play. An example 
in which these phases play a fundamental role was men- 
tioned in the Introduction. In Ref. it was shown that 
during the free expansion of the hard-core boson gas on a 
lattice n k approaches the one of noninteracting fermions. 
This occurs due to the effects of Q x (j) since |Ae(t)| was 
found to decay with exactly the same power law than in 
equilibrium, which would mean that the system should 
have a large peak at k = 0. 




x/a 

FIG. 10: (Color online) Time average of the modulus and 
real part of the one-particle density matrix compared with 
the power-law decay present in the initial state and with the 
results obtained using the usual (thermal) and generalized 
(constrained) Gibbs distributions. We calculated p x from the 
center of the box, and for the time average we measured p x 
in time intervals At — 20t between r = 20t and r = 5000i. 
These systems have 300 lattice sites, 150 particles, and L — 2. 
The initial value of A is A = and the final one A = l.Oi. 
Straight lines following the constrained and thermal results 
signal an exponential decay of p x - 

We find that after relaxation both the modulus and 
the real part of the one-particle density matrix in the 
time evolving system are very different from p x in the 
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thermal and constrained ensemble, for large values of x. 
This can be seen in Fig. [TU] where we have plotted time 
average and thermodynamic results. [We did not include 
results for lm(p x ) since it oscillates strongly between pos- 
itive and negative values.] Figure [TOl shows that while in 
the usual and generalized Gibbs ensemble one-particle 
correlations decay exponentially at large distances (with 
different correlation lengths in each case), no such simple 
exponential decay can be seen cither in the modulus or 
the real part of p x in the evolving system. In addition, 
the decay of the correlations in the latter case cannot 
be described by a power law like the one in the initial 
state. Their decay is actually faster. (We also show in 
the figure p x in the initial superfluid state.) The non- 
trivial effect of the phases x {t) is evident in n&, which 
is the diagonal part of the Fourier transform of the one- 
particle density matrix. While the mean value of \p x \ in 
the out-of-equilibrium system decays more slowly than 
p x in the constrained thermodynamics, the obtained 
from both of them are identical [Fig. [Ha)]. We then 
conclude that there is no simple comparison possible be- 
tween the time average of one-particle correlations in the 
system out of equilibrium (after relaxation) and in its 
"equivalent" generalized Gibb's ensemble. Hence, we do 
not consider that this quantity, which is not a physical 
observable, is an indicator of thermalization. 



IV. HCB'S IN A SUPERLATTICE PLUS A 
HARMONIC CONFINING POTENTIAL 

In this, and the following section, we generalize the re- 
sults obtained in the previous sections to inhomogeneous 
systems. In particular, we study the effects of having 
a harmonic confining potential, relevant to experimental 
systems, in addition to the superlattice potential. In this 
case the Hamiltonian reads 

H = -^(t M+1 6j& i+1 + H.c.) + A^cos( ^ 

i i ^ 

+v 2 ^2^m, (26) 



where V2 is the curvature of the harmonic trap. As in the 
previous sections in what follows we restrict the analysis 
to the case L = 2. 

From earlier studies of the bosonic Hubbard model in 
a trap it is known that in the presence of a confining 
potential superfluid and Mott-insulating phases coexist 
space separated [Ell HH, [53| ■ On the other hand, in the 
HCB case it has been shown that the presence of the trap 
does not destroy the power-law decay of the one-particle 
correlations known from the homogeneous case [38(. 

Previous studies have shown that the key parameter 
that controls the thermodynamic behavior of these con- 
fined systems is the characteristic density 



where £ = (V^/f) is a length scale set by the com- 
bination of the trap and the underlying lattice (38|. As 
p — > 0, one recovers the continuum limit. On the other 
hand, as p increases beyond a critical value insulating 
regions build up in the system. In what follows we nor- 
malize distances in the trap by £ and calculate the mo- 
mentum distribution function as 



n fe = ^£e-^-0,., 



(28) 



In Figs. [TITa)-[TlTd) we show the density profiles (av- 
eraged per unit cell) for harmonically trapped systems 
with A = 0.5i, at different fillings. These plots show 
that like in the Bosc-Hubbard case superfluid and insu- 
lating (constant density n = 0.5 and 1) phases coexist 
space separated in the trap. The width of the insulating 
phases with n — 0.5 is determined by the gap, which for 
the cases depicted in Fig. [11] is A = t. The effect of the 
curvature of the harmonic potential is already considered 
by normalizing the distances by £. 
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FIG. 11: (Color online) Density (averaged per unit cell) (a)- 
(d) and normalized momentum distribution function (e)-(h) 
for trapped systems with 1000 lattice sites, V2CE 2 = 3x 10 _5 t, 
L — 2, and A = 0.5t. In (a)-(d), the regions with constant 
density are local insulators. 

The formation of insulating domains in the trap 
strongly supresses the peaks observed in the momen- 
tum distribution at k = 0. This can be seen in Figs. 
lllf e)- fTTT h). where we have plotted the momentum pro- 
files corresponding to the densities in Figs. [TlTa)-[TTTd'). 
However, these peaks are not destroyed and can be very 
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sharp. As can be seen in Fig. |T2l the reason is that quasi- 
long-range correlations are still present in the superfluid 
regions. We find that their power-law decay is the same 
observed in the absence of the superlattice [38[ and in the 
periodic case discussed in Sec. II — i.e., p x ~ 1/y/x, with 
I . Here of course translational invariance is 
broken by the trap; however, the above power law can be 
seen independently of the points Xi and Xj chosen within 
the superfluid part. On the other hand, as expected, in 
the insulating phases, where the density is constant, the 
decay of p x is exponential, exactly like the one in Fig. 
121 The particles there, being localized, are the ones that 
mainly contribute to the high momentum tails of n*,. 
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FIG. 12: (Color online) Decay of one-particle correlations p x 
(x = |a?i — a?j |) (averaged per unit cell as described in Sec. 
II) in the superfluid domains. We have chosen x% to be the 
center of the trap [see the corresponding density profiles in 
Figs. lllf a') and lllf c')]. Straight lines signal a power-law decay 
px ~ l/\fx, which disappears when rij — > 0, 0.5, or 1. In 
the inset we show the scaling of the lowest natural orbital 
occupation with an increasing number of particles keeping 
constant the characteristic density of the system, which in 
this case is p = 0.66. The straight line shows that Ao ~ V Nb. 

More information about the physics of the trapped sys- 
tem can be obtained studying the natural orbitals ((j> v ), 
which can be considered to be like effective single-particle 
states when the system consists of interacting particles. 
They are defined as the eigenfunctions of the one-particle 
density matrix pij 

N 

\<t>l (29) 

3=1 

A,, denotes the occupation of each orbital. In higher 
dimensions, when only the lowest natural orbital (the 
highest occupied one) scales ~ Nb, it can be regarded 
as the BEC order parameter — i.e., the condensate [55| . 
This scaling of the lowest natural orbital in higher dimen- 
sions is related to the presence of long-range off-diagonal 
order (EH ]. In the one-dimensional case we are study- 
ing here there is no long-range order, only quasi-long- 
range order characterized by a l/y/x decay of the corre- 
lations (Fig. [T2"|) . which as shown in the inset of Fig. [TJ] 
produces a power-law (y/Nb) scaling of the lowest nat- 



ural orbital occupation. (These scaling laws have also 
been observed in harmonically trapped continuous sys- 
tems [H, H^].) Hence, we refer to the lowest nat- 
ural orbital as a quasicondensate since Ao — > oo when 
N b -> oo, but \ /N b -> 0. 
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FIG. 13: (Color online) Natural orbital occupations (a)- 
(d) and normalized wave function of the lower natural or- 
bital ip° = (Nb^/a) 1 ^ 4 (f> (averaged per unit cell) (e)-(h), for 
trapped systems with 1000 lattice sites, V2C1 2 = 3 x 10 -5 t, 
L — 2, and A — 0.5t. In the insets of (a)-(d) we show the 
occupation of the lowest 11 natural orbitals, where degener- 
acy can be seen in the presence of the insulating domains. 
Comparing (e)-(h) with Figs. Illf al- TilT d') one can see that 
the lowest natural orbitals only have a finite weight in the 
superfluid domains. 

In Figs.[T5{a)-{T3Jd) we show the occupation of the nat- 
ural orbitals (ordered from the highest occupied one to 
the lowest occupied one) as a function of the orbital num- 
ber 77. The formation of the insulating domains not only 
reduces the occupation of the lowest natural orbitals, but 
also makes them degenerate (insets). The normalized 
wave function of the lowest natural orbitals [H[ 

/ = (7V fc C/a) 1/4 ^°, (30) 

are plotted in Figs. fl3Te)-|T5yh) . Comparing these wave 
functions with the density profiles in Figs. [TTTa)- [TlT d) 
one can see that these orbitals are localized in the su- 
perfluid regions and have zero weight in the insulating 
domains. Degeneracy then appears because two identical 
quasicondensates develop to the sides of the central insu- 
lating core in the cases depicted in Figs. fTTT b) and[TlTd). 
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According to the number of particles in the system the 
lowest natural orbitals can be located outside or inside 
the insulating domains with average density n = 0.5. 

The behavior of the lowest natural orbital occupations 
with an increasing number of particles in the trap can be 
seen in Fig. [14] It perfectly reflects the formation and 
destruction of insulating domains in the system. The 
first degeneracy for the four lowest natural orbitals (Fig. 
I14[) appear when the insulating phase with mean den- 
sity n = 0.5 sets in the middle of the trap [Fig. ITlT b)]. 
This degeneracy disappears for the natural orbitals 3 and 
4 when a supcrfluid domain, like the one in Fig. Illf c). 
develops in the center of the system. The central su- 
pcrfluid phase grows with increasing number of particles, 
and the quasicondensate there becomes the highest pop- 
ulated. Finally, degeneracy sets up once again when the 
full filled insulator (n = 1) appears in the center of the 
trap [like in Fig. HUd)] ■ 



1 




and A — 0.5t. Degeneracy sets in when insulating domains 
(with density n = 0.5 first and n = 1.0 second) appear in the 
middle of the trap. 



V. DYNAMICS AND RELAXATION IN A 
HARMONIC TRAP 

In this section we study the dynamics in the presence of 
a harmonic trap. Since at very low densities in the trap, 
when the average interparticlc distance is much larger 
than the lattice spacing, the lattice system and the con- 
tinuum are equivalent, we start analyzing this case. 

We perform a numerical experiment similar to that 
recently done at Penn State [171 ] . There, in order to gen- 
erate a highly exited state in a trap, an optical lattice 
potential was applied during a short period of time to 
an array of ID Bose gases. This produced a momen- 
tum distribution with extra peaks away from k = (the 
only one present in the absence of a lattice) . After turn- 
ing off the lattice potential the dynamics of the system 
was studied using timc-of-flight measurements. It was 




1/t 

FIG. 15: (Color online) Evolution of the momentum distri- 
bution function after a superlattice potential with L = 8 and 
A = 2.0t is applied for a short period of time (At = 0.5t) to 
the ground state of trapped system with 10 HCB's, V2<i 2 = 
4 x 10 -6 £, and 2000 lattice sites. Undamped oscillations can 
be seen in during our simulation time. 



found that after some periods of oscillation the system 
relaxed to a momentum distribution that was not the 
one in thermal equilibrium. Relaxation to an equilib- 
rium distribution was attributed mainly to the effects of 
the anharmonicity of the confining potential. In Ref. [lj| 
and in this paper, we have shown why in the integrable 
limit of infinitely strong interactions the system relaxes 
to a distribution that is not the one in thermal equilib- 
rium. Surprisingly, in the experiments (l7| it was found 
that the absence of thermalization extends towards the 
nonintegrable region of finite interactions. Even though 
these ID systems can be very well described by the Lieb- 
Linigcr model [l5[, which is integrable in periodic homo- 
geneous systems, the presence of the trapping potential 
in the experiment breaks integrability but did not allow 
the system to thermalize. 

To recreate the experiment in Rcf. jTJ (performed in 
continuous space) in a perfect harmonic potential, we 
consider a very dilute system in our lattice model. We 
study the dynamics of 10 HCB's in 2000 lattice sites; 
i.e., the interparticle distance is much larger than the 
lattice spacing so that the effects of the underlying lat- 
tice are negligible. We then turn on, for a short period of 
time, an additional periodic potential with L = 8. The 
dynamics of the momentum distribution function after 
turning off this additional lattice is shown in Fig. 1151 
This figure shows that no relaxation towards an equilib- 
rium distribution can be seen in the perfectly harmonic 
case. The additional peaks in nk oscillate back and forth 
almost without damping. This occurs even when many 
eigenmodes of the system have been excited by turning 
on and off the latti ce p otential. Our results confirm the 
conclusion in Rcf. [17|] that relaxation in the hard-core 
regime may be mainly related to the anharmonicity of 
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the confining potential. 

While the low-energy region (low occupation) of the 
harmonic trap spectrum in the presence of the lattice is 
like the one in the continuum, this is not true anymore in 
the high-energy region (high occupation); i.e., the lattice 
model and the continuum one start to differ (38l . l6(| . The 
spectrum in the lattice departs from the linear dispersion 
relation of the harmonic potential, and degeneracies set 
in when insulating domains develop in the ground state 
[6lj . Hence, for large fillings in the trap we may expect 
to see relaxation towards an equilibrium distribution. 

We study the dynamics of systems in which initially, in 
the presence of a superlattice potential, there is a coexis- 
tence of superfiuid and insulating domains. We then turn 
off the superlattice potential and let the system evolve in 
the presence of the trap. In Fig. [I])] we show the evolution 
of the occupations of the zero- momentum state (nk=o) 
and the lowest natural orbital (Ao) when (a) the initial 
state has a half-filled insulator in the center of the trap 
[Fig. UTTa)] and (b) two insulating shoulders surround a 
central superfiuid region [Fig. [TTT d)]. Figure [TrJl shows 
that for high fillings there is a strong damping of the os- 
cillations of rik=o an d Ao, in contrast to the harmonically 
trapped case without the lattice. Hence, in experiments 
relaxation to an equilibrium state is to be expected when 
a lattice is present along the ID tubes. 
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FIG. 16: (Color online) Evolution of the occupation of the 
zero-momentum state (top plots) and the occupation of the 
lowest natural orbital (bottom plots) after a superlattice po- 
tential is turned off in trapped systems with 900 lattice sites 
and V^a 2 = 3 x 10 -5 t. The evolution starts from the ground 
state of a system with L = 2 and A = 0.5t. The number of 
particles is N b = 200 (a) and Ni = 299 (b). The correspond- 
ing initial density profiles can be seen in Figs.[T7Ja) and !17f dl. 
The dashed lines are the results obtained with the generalized 
Gibbs distribution explained in the text. 

As seen in Fig. \W\ after t = 5000i the oscillations of 
rik=o and Ao are very small and occur around an approx- 
imately constant value. In Fig. 1171 we compare the time 



average (between r = 5000 and lOOOOi) of the density 
profiles, the momentum distribution function, and the 
natural orbitals, with the predictions of the usual grand- 
canonical ensemble and the generalized Gibbs distribu- 
tion introduced in Ref. [ll| . For all quantities the results 
of the time averages and the generalized Gibbs distribu- 
tion are on top of each other. On the other hand, the 
differences between time averages and the thermal en- 
semble are apparent in all cases and particularly large in 
the low-momentum region of and in the occupation 
of the lowest natural orbitals. Only in the density pro- 
files are the differences smaller and mainly visible in the 
regions where the density approaches zero. 

Since the density profiles of HCB's and noninteracting 
fcrmions are identical and their time average coincides 
with the results of the generalized Gibbs distribution, 
one realizes that noninteracting systems are the simplest 
case to which the constrained thermal equilibrium theory 
[T§ | explained in Sec. Ill can be applied. 
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FIG. 17: (Color online) Time average of the density profiles 
(a),(d), momentum distribution function (b),(e), and the oc- 
cupation of the lowest 400 natural orbitals. The average is 
performed between r — 5000i and r = lOOOOi (see Fig. [T6)l 
with measurements done in time intervals At = 40t. The 
evolution is performed in trapped systems with 900 lattice 
sites, V2CL 2 — 3 x W~ 5 t, starting from the ground state in the 
presence of a superlattice with L = 2 and A = 0.5t, after 
turning off the superlattice. The results of the time average 
are compared with the ones obtained in the usual thermal en- 
semble and the constrained theory explained in the text. The 
number of particles is N b = 200 (a)-(c) and N b = 299 (d)- 
(f). In (a) and (d) we included the averaged density per unit 
cell in the initial state. Flat regions correspond to insulating 
domains. 
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VI. CONCLUSIONS 

We have studied ground-state properties and the 
noncquilibrium dynamics of hard-core bosons on one- 
dimensional lattices in the presence of an additional pe- 
riodic potential (supcrlattice) and a harmonic trap. 

In the periodic case the superlattice potential opens 
gaps at the borders of the reduced Brillouin zones gener- 
ating insulating phases with fractional fillings. In these 
insulating phases one-particle correlations decay expo- 
nentially [p x ~ exp(— x/£)], and we have studied how the 
correlation length £ behaves as a function of the strength 
A of the supcrlattice potential. On the other hand, we 
find that in the gapless superfluid phases one-particle 
correlations exhibit quasi-long-range order. They decay 
with the same power law shown in Ref. [38| to be univer- 
sal in the absence of the superlattice potential. 

In the presence of an additional confining potential, 
which we have taken to be harmonic for a closer connec- 
tion to the experiments, superfluid and insulating do- 
mains coexist phase separated. We have shown that 
in the superfluid domains, where the density changes 
spatially, one-particle correlations decay with the same 
power law than in the homogeneous periodic case. This 
decay produces a square-root scaling of the occupation 
of the lowest natural orbital with the number of particles 
in the trap and independently of the presence or not of 
insulating domains. 

We have also studied the noncquilibrium dynamics of 
these systems after a quench of the supcrlattice poten- 
tial A. In particular, we have considered a half-filled box 
in the presence of a superlattice with period 2, since for 
A = the system is superfluid while for A ^ it is insu- 
lating. After a sudden switch-on of A we have shown that 
the initial momentum distribution (rife) of the superfluid 
phase collapses and revives with a period determined by 
A, like in the experiment in Ref. [501 ] where the period 
was determined by the on-site repulsion U. After several 
oscillations, the number depending on A and the hopping 
parameter t, we have also seen that the system relaxes 
to an equilibrium distribution with very small fluctua- 
tions in rife. The time average of this physical observable 



was then shown to be very well described by a general- 
ized Gibbs distribution introduced in Ref. [l8|. On the 
other hand, after a sudden switch-off of A we have found 
that not only the generalized Gibbs ensemble but also 
the usual thermal ensemble describes very well the time 
average result when the system is initially deep in the 
insulating regime (A > i). This is indeed a surprising 
result since we are seeing thcrmalization in an integrable 
system. But of course it is only seen in the very particu- 
lar case in which the quench of the interaction drives the 
system from an insulating state to a superfluid one. In 
the latter the temperature plays a very similar role than 
the gap present in the initial insulator. In transitions 
between different superfluid states [l8[ no such thermal- 
ization is seen. 

Finally, we have considered the nonequilibrium dynam- 
ics for the harmonically confined case. We have shown 
that in this more experimentally relevant system the time 
average of obscrvablcs, like density and momentum dis- 
tribution, can also be very well described by the gener- 
alized Gibbs distribution. Damping occurs in this case, 
even in a perfect harmonic trap, because of the presence 
of the lattice. At very low densities, equivalent to the 
continuum case, we have shown that oscillations of 
occur almost without damping. 
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